import numpy as np

erf_table = np.array([   
           0,   578,  1155,  1732,  2308,  2883,  3456,  4028,  4598,
        5165,  5731,  6293,  6852,  7408,  7961,  8510,  9055,  9595,
       10131, 10663, 11189, 11711, 12227, 12737, 13242, 13741, 14234,
       14721, 15201, 15675, 16142, 16602, 17056, 17502, 17941, 18373,
       18798, 19215, 19625, 20028, 20422, 20809, 21189, 21561, 21925,
       22281, 22629, 22970, 23303, 23628, 23946, 24256, 24558, 24853,
       25140, 25420, 25693, 25958, 26215, 26466, 26709, 26946, 27175,
       27398, 27614, 27823, 28026, 28222, 28412, 28596, 28773, 28945,
       29111, 29271, 29425, 29574, 29718, 29856, 29990, 30118, 30242,
       30360, 30474, 30584, 30689, 30791, 30888, 30981, 31070, 31155,
       31237, 31315, 31390, 31461, 31530, 31595, 31657, 31717, 31774,
       31828, 31879, 31928, 31975, 32019, 32062, 32102, 32140, 32176,
       32211, 32243, 32274, 32303, 32331, 32358, 32382, 32406, 32428,
       32449, 32469, 32488, 32506, 32522, 32538, 32553, 32567, 32580,
       32592, 32604, 32615, 32625, 32635, 32644, 32652, 32660, 32668,
       32674, 32681, 32687, 32693, 32698, 32703, 32708, 32712, 32716,
       32720, 32724, 32727, 32730, 32733, 32735, 32738, 32740, 32742,
       32744, 32746, 32748, 32749, 32751, 32752, 32754, 32755, 32756,
       32757, 32758, 32758, 32759, 32760, 32761, 32761, 32762, 32762,
       32763, 32763, 32764, 32764, 32764, 32765, 32765, 32765, 32765,
       32766, 32766, 32766, 32766, 32766, 32767, 32767, 32767, 32767,
       32767, 32767, 32767, 32767, 32767, 32767, 32767, 32768, 32768,
       32768, 32768, 32768, 32768, 32768, 32768, 32768, 32768, 32768,
       32768, 32768, 32768, 32768, 32768, 32768, 32768, 32768, 32768,
       32768, 32768, 32768, 32768, 32768, 32768, 32768, 32768, 32768,
       32768, 32768, 32768, 32768, 32768, 32768, 32768, 32768, 32768,
       32768, 32768, 32768, 32768, 32768, 32768, 32768, 32768, 32768,
       32768, 32768, 32768, 32768, 32768, 32768, 32768, 32768, 32768,
       32768, 32768, 32768, 32768], dtype=np.int32)

def np_erf(x):
    # save the sign of x
    sign = ((x >= 0).astype(int) * 2) - 1
    x = np.abs(x)

    # constants
    a1 =  0.254829592
    a2 = -0.284496736
    a3 =  1.421413741
    a4 = -1.453152027
    a5 =  1.061405429
    p  =  0.3275911

    # A&S formula 7.1.26
    t = 1.0/(1.0 + p * x)
    y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*np.exp(-x*x)
    return sign*y # erf(-x) = -erf(x)

def erf_lut(x, nearest=True):
    abs_x = np.abs(x.astype(np.int32)) >> 7  # input in Q13
    abs_x_masked = np.where(abs_x >= 255, 0, abs_x)
    if nearest:
        result = np.where(abs_x >= 255,
                          0x7FFF,
                          erf_table[abs_x_masked])
    else:
        ua = erf_table[abs_x_masked]
        ub = erf_table[abs_x_masked+1]
        ut = abs_x & 0xFF
        result = np.where(abs_x >= 255,
                          0x7FFF,
                          ((ua << 8) + ut * (ub-ua)) >> 8) # Q15*Q8 = Q24
    return np.where(x > 0, result, -result)
